/********************************************************************************

    C++ Implementation of Symmetric Level-Index Arithmetic (see readme.txt for usage)
    Copyright (C) 2007  Xunyang Shen
    Email: xunyang.shen@gmail.com

    This library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
    License as published by the Free Software Foundation; either
    version 2.1 of the License, or (at your option) any later version.

    This library is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with this library; if not, write to the Free Software
    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

**********************************************************************************/
/*
    If you are going to make contributions to the code, keep in mind that
    performance is the first priority.  The "bad" OO design in this library
    is not a defect here.  On the contrary, long function bodies make the
    implementation runs effiently, because overhead is avoided as much as
    possible.

**********************************************************************************/

#include "SLI.h"

#define TRUE 1
#define FALSE 0
#define INTERVAL 32
#define INTLENGTH 0.03125
#define FLPMAX 6.7039039649712985e+153  // 2^511
#define FLPMIN 1.4916681462400413e-154  // 2^-511
#define SLIMIN 4.57088361887805  // psi(2^511)
#define REALMAX 1.7976931348623157e+308
#define REALMIN 2.2250738585072014e-308
#define HYPERMAX 7.104758600739439e+002
#define FLOATINF 1e400
#define HSEQUENCE while (h >= 1) {h = log(h); j++;}
#define ABS(x) ((x) >= 0 ? (x) : -(x))

using namespace std;

//////////////////////////////////////////////////////////////////////
// Constructors
//////////////////////////////////////////////////////////////////////

inline SLI::SLI(unsigned int slevel, double sindex, unsigned int ssignl, unsigned int ssigni)
{
	level = slevel;
	index = sindex;
	signl = ssignl;
	signi = ssigni;

	if (slevel > 7)	// overflow or underflow
	{
		level = 7;
		index = 1; 
	}

	if (slevel + sindex < SLIMIN) // in between 2^-511 and 2^511
	{
		index = (*this).Double();
		level = 0;
		signl = 0;
		signi = 0;
	}
}

SLI::SLI(short slevel, double sindex) throw(domain_error)
{
	if (ABS(slevel) > 7)
	{
		cout << "The level must be an integer between -7 and 7!" << endl;
		system("PAUSE");
		throw domain_error("Level domain violation");
	}

	if (slevel >= 0)
	{
		level = slevel;
		signl = 0;
	}
	else
	{
		level = -slevel;
		signl = 1;
	}

	if (sindex >= 0)
	{
		index = sindex;
		signi = 0;
	}
	else
	{
		index = -sindex;
		signi = 1;
	}

	if (index >= 1)
	{
		cout << "The index must be a number between -1 and 1!" << endl;
		system("PAUSE");
		throw domain_error("Index domain violation");
	}

	if (level + index < SLIMIN) // in between 2^-511 and 2^511
	{
		index = (*this).Double();
		level = 0;
		signl = 0;
		signi = 0;
	}
}

SLI::SLI(double a)
{
	if (a >= 0)
	{
		index = a;
		signi = 0;
	}
	else
	{
		index = -a;
		signi = 1;
	}

	if ((index > FLPMIN) && (index < FLPMAX))  // use FLP
	   //pow(2, -511) < abs(a) < pow(2, 511)
	   //511 * 2 = 1022 < 1024, which guarentees pure FLP multiplications work properly
	{
		level = 0;
		index = a;
		signl = 0;
		signi = 0;
	}

	else if (a == 0)  // 0
	{
		level = 7;
		index = 1;
		signl = 1;
		signi = 0;
	}

	else if (index <= REALMAX) // use SLI
	{
		level = 0;
		if (index < REALMIN)  //gradual underflow in FLP
			index = REALMIN;

		if (index >= 1) // large
			signl = 0;
		else   //small
		{
			signl = 1;
			index = -log(index);
			level++;
		}
	
		while (index >= 1)
		{
			index = log(index);
			level++;
		}
	}

	else
	{
		level = 7;
		index = 1;
		signl = 0;
	}
}

SLI::SLI(double a, char b)
{
	if (a >= 0)
	{
		index = a;
		signi = 0;
	}
	else
	{
		index = -a;
		signi = 1;
	}

	level = 0;
	
	if (index >= 1) // large SLI
		signl = 0;
	else   //small SLI
	{
		signl = 1;
		index = -log(index);
		level++;
	}
	
	while (index >= 1)
	{
		index = log(index);
		level++;
	}
}

SLI::SLI()  // initialize as 0
{
	level = 7;
	index = 1;
	signl = 1;
	signi = 0;
}

//////////////////////////////////////////////////////////////////////
// Turn a SLI number to a decimal floating-point number
// and return the exponent or significand
//////////////////////////////////////////////////////////////////////
double SLI::GetDecimalExp() const
{
	if (level == 0)  //FLP
		return floor(log10(ABS(index)));

	// the exponent is too big or too small to be represented in double precision FLP arithmetic
	if ((level > 5) || ((level == 5) && (index > 0.6322952489608)))
		if (signl)
			return 0;
        else 
			return FLOATINF;

	double log10_X = (SLI(level-1, index, 0, signl) / log(10.0)).Double();
	return floor(log10_X);
}

double SLI::GetDecimalSig() const
{
	if (level == 0)  //FLP
		return pow(10.0, fmod(log10(index), 1));

	// the exponent is too big or too small to be represented in double precision FLP arithmetic
	if ((level > 5) || ((level == 5) && (index > 0.6322952489608)))
		if (signl)
			return 0;
        else 
			return signi ? -FLOATINF : FLOATINF;
			
	double log10_X = (SLI(level-1, index, 0, signl) / log(10.0)).Double();
	double f_X = pow(10.0, fmod(log10_X,1));
	return signi ? -f_X : f_X;	
}

//////////////////////////////////////////////////////////////////////
// Displaying a SLI number as
// (+/-level, +/-index)
//////////////////////////////////////////////////////////////////////

void display(const SLI& s)
{
	cout << setprecision(15);
    if (s.level == 0)  //FLP
		cout << s.index;
	else  //SLI
	{
		if ((s.index == 1) && (s.level == 7)) //overflow or underflow
		{
			if (s.signl)
				cout << 0 * 1e-10;
			else if (s.signi)
				cout << "-Inf";
			else
				cout << "Inf";
		}
		else
		{
			cout << "(" << (s.signl ? "-" : "") << s.level << ", ";
            cout << (s.signi ? "-" : "") << s.index << ")";
        }
	}
}

//////////////////////////////////////////////////////////////////////
// Display a SLI number in the form of
// Mantissa * 10 ^ Exponent
//////////////////////////////////////////////////////////////////////
void dispfloat(const SLI& s)
{
	cout << setprecision(15);
	if (s.level == 0)  //FLP
		cout << s.index;
	else  //SLI
		// the exponent is too big or too small to be represented in double precision FLP arithmetic
		if ((s.level > 5) || ((s.level == 5) && (s.index > 0.6322952489608)))
		{
			if (s.signl)
				cout << 0 * 1e-10;
			else if (s.signi)
				cout << "-Inf";
			else
				cout << "Inf";
		}
		else
		{
			double log10_X = (SLI(s.level-1, s.index, 0, s.signl) / log(10.0)).Double();
			double e_X = floor(log10_X);
			double f_X = pow(10.0, double(log10_X - e_X));

			if (s.signi)
				f_X = -f_X;
			cout << f_X << " * 10 ^ " << e_X;
		}
}

//////////////////////////////////////////////////////////////////////
// Get the level/index of a SLI number
//////////////////////////////////////////////////////////////////////
unsigned SLI::GetLevel() const
{
	if (level == 0)
		return SLI(index, 'S').level;
	else
		return level;
}

unsigned SLI::GetSignl() const
{
	if (level == 0)
		return SLI(index, 'S').signl;
	else
		return signl;
}

unsigned SLI::GetSigni() const
{
	if (level == 0)
		return SLI(index, 'S').signi;
	else
		return signi;
}

double SLI::GetIndex() const
{
	if (level == 0)
		return SLI(index, 'S').index;
	else
		return index;
}

//////////////////////////////////////////////////////////////////////
// Get the signed level/index
//////////////////////////////////////////////////////////////////////
short SLI::GetSignedLevel() const
{
	if (level == 0)
	{
		SLI s(index, 'S');
		return (s.signl ? -short(s.level) : s.level);
	}
	else
		return (signl ? -short(level) : level);
}

double SLI::GetSignedIndex() const
{
	if (level == 0)
	{
		SLI s(index, 'S');
		return (s.signi ? -s.index : s.index);
	}
	else
		return (signi ? -index : index);
}

//////////////////////////////////////////////////////////////////////
// Convert to double
//////////////////////////////////////////////////////////////////////
double SLI::Double() const
{
	double res;
	unsigned i;
	
	if ((index == 1) && (level == 7) && signl)
	   return 0;
	
	res = index;
	for (i=1; i<=level; i++)
		res = exp(res);
	res = (signl ? 1/res : res);
	return (signi ? -res : res);
}

//////////////////////////////////////////////////////////////////////
// Absolute value
//////////////////////////////////////////////////////////////////////
SLI abs(const SLI& s)
{
	if (s.level == 0)
		return s.index >= 0 ? SLI(s.index) : SLI(-s.index);
	else
		return SLI(s.level, s.index, s.signl, 0);
}

//////////////////////////////////////////////////////////////////////
// Reciprocal
//////////////////////////////////////////////////////////////////////
SLI SLI::Rec() const
{
	if (level == 0)
		return SLI(1/index);
	else
		return SLI(level, index, ~signl, signi);
}

//////////////////////////////////////////////////////////////////////
// Unary minus
/////////////////////////////////////////////////////////////////////
SLI SLI::operator -() const
{
	if (level == 0)
		return SLI(-index);
	else
		return SLI (level, index, signl, ~signi);
}

//////////////////////////////////////////////////////////////////////
// Exponential function
//////////////////////////////////////////////////////////////////////
SLI exp(const SLI& s)
{
	if (s.level == 0)
		if (ABS(s.index) <= 708) // prevent overflow
			return SLI(exp(s.index));
		else
			return exp(SLI(s.index, 'S'));
	else if (s.signl)
		return SLI(1.0);
	else
		return SLI(s.level+1, s.index, s.signi, 0);
	       //increase level, ...  , reciprocal decided by sign, always positive
}

//////////////////////////////////////////////////////////////////////
// Logarithm functions
//////////////////////////////////////////////////////////////////////
SLI log(const SLI& s)
{
	if (((s.level != 0) && (!s.signi)) || ((s.level == 0) && (s.index > 0))) // None negative
	{
		if (s.level == 0)
			return SLI(log(s.index));
		else if (s.index == 1) // 0 or inf
			return SLI(7, 1, 0, s.signl);
		else
			return SLI(s.level-1, s.index, 0, s.signl);
				//decrease level, ..., 0 for level beyond 1/-1, sign decided by reciprocal  
	}
	else
		return SLI(7, 1, 0, 1); // -Inf
}

SLI log2(const SLI& s)
{
	static double LN2 = log(2.0);
	if (((s.level != 0) && (!s.signi)) || ((s.level == 0) && (s.index > 0))) //none negative
	{
		if (s.level == 0)
			return SLI(log(s.index) / LN2);
		else if (s.index == 1) // 0 or inf
			return SLI(7, 1, 0, s.signl);
		else
			return SLI(s.level-1, s.index, 0, s.signl) / LN2;
				//decrease level, ..., 0 for level beyond 1/-1, sign decided by reciprocal  
	}
	else
		return SLI(7, 1, 0, 1); // -Inf
}

SLI log10(const SLI& s)
{
	static double LN10 = log(10.0);
	if (((s.level != 0) && (!s.signi)) || ((s.level == 0) && (s.index > 0))) //none negative
	{
		if (s.level == 0)
			return SLI(log10(s.index));
		else if (s.index == 1) // 0 or inf
			return SLI(7, 1, 0, s.signl);
		else
			return SLI(s.level-1, s.index, 0, s.signl) / LN10;
				//decrease level, ..., 0 for level beyond 1/-1, sign decided by reciprocal  
	}
	else
		return SLI(7, 1, 0, 1); // -Inf
}

//////////////////////////////////////////////////////////////////////
// Relational operators
// In order to reduce the overhead of function calls, some code are
// copied in different functions instead of using subroutines.
// So make sure to CHANGE ALL OF THEM if any modification is made!
//////////////////////////////////////////////////////////////////////
bool operator ==(const SLI& p, const SLI& q)
{
	return ((p.level == q.level) && (p.index == q.index) && (p.signl == q.signl) && (p.signi == q.signi));
}

bool operator !=(const SLI& p, const SLI& q)
{
	return ((p.level != q.level) || (p.index != q.index) || (p.signl != q.signl) || (p.signi != q.signi));
}

bool operator >(const SLI& p, const SLI& q)
{
	if (p.level == 0)
		if (q.level == 0)
			return (p.index > q.index);
		else
			return (SLI(p.index, 'S') > q);
	if (q.level == 0)
		return (p > SLI(q.index, 'S'));
		
	if (p.signi != q.signi)  //different sign
		return (p.signi < q.signi);

	bool res;
	short lv1, lv2;
	lv1 = (p.signl ? -short(p.level) : p.level);
	lv2 = (q.signl ? -short(q.level) : q.level);

	if (lv1 == lv2)  //same level
		res = ((p.signl ? -p.index : p.index) > (q.signl ? -q.index : q.index));
	                 //compare indices combining signs of level
	else  //different leve
		res = (lv1 > lv2);

	return (p.signi ? !res : res);  //return uminus for a negative number
}

bool operator <(const SLI& p, const SLI& q)
{
	if (p.level == 0)
		if (q.level == 0)
			return (p.index < q.index);
		else
			return (SLI(p.index, 'S') < q);
	if (q.level == 0)
		return (p < SLI(q.index, 'S'));

	if (p.signi != q.signi)  //different sign
		return (p.signi > q.signi);

	bool res;
	short lv1, lv2;
	lv1 = (p.signl ? -short(p.level) : p.level);
	lv2 = (q.signl ? -short(q.level) : q.level);

	if (lv1 == lv2)  //same level
		res = ((p.signl ? -p.index : p.index) < (q.signl ? -q.index : q.index));
	                 //compare indices combining signs of level
	else  //different leve
		res = (lv1 < lv2);

	return (p.signi ? !res : res);  //return uminus for a negative number
}

bool operator >=(const SLI& p, const SLI& q)
{
	if (p.level == 0)
		if (q.level == 0)
			return (p.index >= q.index);
		else
			return (SLI(p.index, 'S') >= q);
	if (q.level == 0)
		return (p >= SLI(q.index, 'S'));

	if (p.signi != q.signi)  //different sign
		return (p.signi < q.signi ? TRUE : FALSE);

	bool res;
	short lv1, lv2;
	lv1 = (p.signl ? -short(p.level) : p.level);
	lv2 = (q.signl ? -short(q.level) : q.level);

	if (lv1 == lv2)  //same level
		res = ((p.signl ? -p.index : p.index) >= (q.signl ? -q.index : q.index));
	                 //compare indices combining signs of level
	else  //different leve
		res = (lv1 > lv2);

	return (p.signi ? !res : res);  //return uminus for a negative number
}

bool operator <=(const SLI& p, const SLI& q)
{
	if (p.level == 0)
		if (q.level == 0)
			return (p.index <= q.index);
		else
			return (SLI(p.index, 'S') <= q);
	if (q.level == 0)
		return (p <= SLI(q.index, 'S'));

	if (p.signi != q.signi)  //different sign
		return (p.signi > q.signi);

	bool res;
	short lv1, lv2;
	lv1 = (p.signl ? -short(p.level) : p.level);
	lv2 = (q.signl ? -short(q.level) : q.level);

	if (lv1 == lv2)  //same level
		res = ((p.signl ? -p.index : p.index) <= (q.signl ? -q.index : q.index));
	                 //compare indices combining signs of level
	else  //different leve
		res = (lv1 < lv2);

	return (p.signi ? !res : res);  //return uminus for a negative number
}

//////////////////////////////////////////////////////////////////////
// Arithmetic operators
//////////////////////////////////////////////////////////////////////
SLI operator +(const SLI& p, const SLI& q)
{
	//Divide SLI interval [4.5625, 5.21875] to 20 subintervals, each subinterval has a length of 1/32.
	//The reason to choose 4.5625 is 4.5625 < 4.57088, which is 2^511
	//The reason to choose 5.21875 is 5.21875 > 5.21574, which makes all additions trivial
	//FL[] records the largest values of y for which 1st order Taylor approximation gives accurate results,
	//corresponding to the x values of the start points of subintervals.
	static double const FL[21] =
		{4.5582433623580281, 4.5909633267988159, 4.6232450529862588, 4.6552150795679417, 4.6869049978709176,
		 4.7184196838818462, 4.7498292887475095, 4.7811659045113117, 4.8124610753762145, 4.8437328906498998,
		 4.8749933766908473, 4.9062475940345553, 4.9374992019728552, 4.9687497726405603, 4.9999999468243512,
		 5.0312499850982526, 5.0624999985163077, 5.0937499998253406, 5.1249999999887086, 5.1562499999995008,
		 5.1874999999999920};
	//DL[] gives the slopes of each interval for linear interpolation.
	static double const DL[21] =
		{1.0470388621052109, 1.0330152379981712, 1.0230408506138531, 1.0140773856952308, 1.0084699523497136,
		 1.0051073557012273, 1.0027717044416704, 1.0014454676768878, 1.0006980887579289, 1.0003355533103218,
		 1.0001349549986571, 1.0000514540255949, 1.0000182613665629, 1.0000055738813103, 1.0000012247648442,
		 1.0000004293777636, 1.0000000418890522, 1.0000000052277755, 1.0000000003453522, 1.0000000000157172,
		 1.0000000000002558};
	//Do similar things to SLI interval [-5.21875, -4.5625]
	static double const FS[21] = 
		{-4.5666055711971163, -4.5963998559046706, -4.6266965867909953, -4.6572728165884403, -4.6880898585426589,
		 -4.7190739395659840, -4.7501717909201169, -4.7813350544478812, -4.8125393838410009, -4.8437669636228007,
		 -4.8750067061452755, -4.9062524146875859, -4.9375007713097911, -4.9687502178398910, -5.0000000540313385,
		 -5.0312500102792708, -5.0625000015194068, -5.0937500001648512, -5.1250000000115410, -5.1562500000004325,
		 -5.1875000000000044};
	static double const DS[21] = 
		{-0.95341711064173751, -0.96949538836238958, -0.97843935351824030, -0.98614534253499642, -0.99149059274640194,
		 -0.99513124333225278, -0.99722443288845852, -0.99853854057982971, -0.99928255301759350, -0.99967176071919539,
		 -0.99986267335393109, -0.99994741191056846, -0.99998228896319574, -0.99999475812631999, -0.99999859993383211,
		 -0.99999971968435375, -0.99999995665422148, -0.99999999509407189, -0.99999999964452968, -0.99999999998630074,
		 -0.99999999999985789};

	if (p.level == 0)  // FLP number
		if (q.level == 0) // pure FLP addition
			return SLI(p.index + q.index);
		else  // mixed SLI and FLP
			return q + p;
	
	if (q.level == 0)  //mixed SLI and FLP
	{
		if (p.signl)
			return SLI(p.Double() + q.index);
	}
	else
	{
		if ((p.signl > q.signl) || ((p.signl == q.signl) && (((!p.signl) && ((p.level < q.level) || ((p.level == q.level) && (p.index < q.index)))) ||
			                                             ((p.signl) && ((p.level > q.level) || ((p.level == q.level) && (p.index > q.index)))))))
			//compare absolute values
			return q + p;
		if ((p.index == q.index) && (p.level == q.level) && (p.signl == q.signl) && (p.signi != q.signi))  // *this == -s
			if ((p.index == 1) && !p.signl)  // +/-Inf -/+ Inf
				return SLI(7, 1, 0, p.signi & q.signi);  // +/- Inf
				// Inf + Inf = Inf, Inf - Inf = Inf, -Inf + Inf = Inf, -Inf - Inf = -Inf
			else
				return SLI(7, 1, 1, 0); // 0
	}

	if ((p.level > 5) || ((p.level == 5) && (p.index > 0.21574365941012072)))  //must be trivial addition
		return p;

	unsigned qsigni, qsignl;
	// index sign of q, level sign of q
	if (q.level == 0)  //mixed SLI and FLP
	{
		qsigni = (q.index < 0);
		qsignl = (ABS(q.index) < 1);
	}
	else
	{
		qsigni = q.signi;
		qsignl = q.signl;
	}

	register unsigned j;
	double a[4], b, c, d, h;

	// Generate sequence a
	a[p.level-2] = exp(-p.index);
	for (j=p.level-2; j>0; j--)
		a[j-1] = exp(-1/a[j]);

	//Calculate d
	if (q.level == 0)
		//if (p.signl == 0) // large or mixed
		d = (q.index >= 0 ? q.index : -q.index) * exp(-1/a[0]);
		//There is no small case in mixed SLI & FLP
		//	d = abs(q.index) * exp(1/a[0]);
	else
	{
		b = exp(-q.index);
		for (j=q.level-2; j>0; j--)
			b = exp(-1/b);
		if (!p.signl)
			if (!qsignl) //large
				d = exp((a[0]-b)/(a[0]*b));
			else //mixed
				d = exp((-a[0]-b)/(a[0]*b));
		else //small
			d = exp((b-a[0])/(a[0]*b));
	}
	
	bool UseTaylor = 0;
	int interNo;
	double intpart, fraction, crit;

	if (p.level + p.index >= 4.5625) 
		if (p.signl)
		{
			fraction = modf((p.level + p.index - 4.5625) * INTERVAL, &intpart) * INTLENGTH;
			interNo = int(intpart);
			crit = FS[interNo] + fraction * DS[interNo];
			if (q.level + q.index > -crit)
				UseTaylor = 1;
		}
		else
			if (q.level == 0)
			{
				if (p.index > 0.5744) // x value could use Taylor for all y < 2^511
					UseTaylor = 1;
			}
			else
			{
				fraction = modf((p.level + p.index - 4.5625) * INTERVAL, &intpart) * INTLENGTH;
				interNo = int(intpart);
				crit = FL[interNo] + fraction * DL[interNo];
				if ((q.signl) || ((!q.signl) && (q.level + q.index < crit)))
					UseTaylor = 1;
			}

	if (UseTaylor)  // 1st order Taylor approximation
	{
		double epsi;
		if ((qsignl) && (p.signl))
			epsi = -d / (d+1); // small
		else
			epsi = d;  // large or mixed
		for (j=0; j<p.level-1; j++)
			epsi *= a[j];

		if (qsigni != p.signi)  // subtraction
			epsi = -epsi;

		double t = p.index + epsi;
		if (t < 0)
			return SLI(p.level-1, exp(t), p.signl, p.signi);
		else if (t > 1)
			return SLI(p.level+1, log(t), p.signl, p.signi);
		else
			return SLI(p.level, t, p.signl, p.signi);
	}
	else  // Full algorithm
	{	// Initial value of sequence c is 1+d, here use log to save storage
		if (p.signi == qsigni)
			c = log(1 + d);
		else
			c = log(ABS(1 - d));
		// Theoretically this ABS() should be removed.
		// But when optimization is used in compiling, 1 - exp((a[0]-b)/(a[0]*b)) could be -0
		// This must be a bug in C++ compiler (both in Dev-C++ and Visual C++)

		// Generate sequence c
		c = p.signl ? 1 - a[0] * c : 1 + a[0] * c;
		j = 0;
		while ((c >= a[j]) && (j < p.level-2))
		       c = 1 + a[++j] * log(c);
		if (c < a[j])
			return SLI(j+1, c/a[j], p.signl, p.signi);
		j++;
		h = p.index + log(c);
		HSEQUENCE
		return SLI(++j, h, p.signl, p.signi);
	}
}

SLI operator -(const SLI& p, const SLI& q)
{
	if (q.level == 0)
		return (p + SLI(-q.index));
	else
		return (p + SLI(q.level, q.index, q.signl, ~q.signi));
}

SLI operator *(const SLI& p, const SLI& q)
{
	if (p.level == 0)  // FLP number
		if (q.level == 0) // pure FLP multiplication
			return SLI(p.index * q.index);
		else  // mixed SLI and FLP
			return q * p;
	
	if ((p.index == q.index) && (p.level == q.level) && (p.signl != q.signl)) // *this == +/- 1/s
		if (p.index == 1) // 0 * (+/-)Inf
			return SLI(7, 1, 0, p.signi ^ q.signi); // +/- Inf
		else
			return SLI(1, 0, 0, p.signi ^ q.signi); // +/- 1

	if ((p.index == 1) && (p.signl)) // 0
		return p;
	
	if ((p.level < q.level) || ((p.level == q.level) && (p.index < q.index)))
		//compare values without signs
		return q * p;

	if ((p.level > 6) || ((p.level == 6) && (p.index > 0.21574365941012072)))  //must be trivial multiplication
		if ((q.index < 0) || q.signi)
			return -p;
		else
			return p;

	unsigned sign;
	register unsigned j;
	double a[4], b, c, c1, h;

	// Generate sequence a
	a[p.level-3] = exp(-p.index);
	for (j=p.level-3; j>0; j--)
		a[j-1] = exp(-1/a[j]);

	// Initial c-sequence c1
	if (q.level == 0)
	{
		sign = p.signi ^ (q.index < 0);
		c1 = p.signl ? 1 - exp(-1/a[0]) * log(ABS(q.index)) : 1 + exp(-1/a[0]) * log(ABS(q.index));
	}
	else
	{
		sign = p.signi ^ q.signi;
		b = exp(-q.index);
		for (j=q.level-3; j>0; j--)
			b = exp(-1/b);
		if (p.signl == q.signl)
			c1 = 1 + exp((a[0]-b)/(a[0]*b));
		else
			c1 = ABS(1 - exp((a[0]-b)/(a[0]*b)));
		// Theoretically this ABS() should be removed.
		// But when optimization is used in compiling, 1 - exp((a[0]-b)/(a[0]*b)) could be -0
		// This must be a bug in C++ compiler (both in Dev-C++ and Visual C++)
	}

	c = 1 + a[0] * log(c1);
	j = 0;
	while ((c >= a[j]) && (j < p.level-3))
		c = 1 + a[++j] * log(c);

	if (c < a[j])
	{
		if (j == 0)
		{
			double a1 = exp(-1/a[0]);
			if (c1 < a1)
				return SLI(1, c1/a1, p.signl, sign);
		}
		return SLI(j+2, c/a[j], p.signl, sign);
	}
	h = p.index + log(c);
	HSEQUENCE
	return SLI(j+3, h, p.signl, sign);
}

SLI operator /(const SLI& p, const SLI& q)
{
	if ((p.level==0) && (p.index==1)) //reciprocal
		return q.Rec();
	if ((p.level == 0) && (q.level == 0))
		return SLI(p.index / q.index);

	if (q.level == 0)
		return (p * SLI(1/q.index));
	else
		return (p * SLI(q.level, q.index, ~q.signl, q.signi));
}

SLI& SLI::operator +=(const SLI& s)
{
	*this = *this + s;
	return *this;
}

SLI& SLI::operator -=(const SLI& s)
{
	*this = *this - s;
	return *this;
}

SLI& SLI::operator *=(const SLI& s)
{
	*this = *this * s;
	return *this;
}

SLI& SLI::operator /=(const SLI& s)
{
	*this = *this / s;
	return *this;
}

SLI pow(const SLI& p, const SLI& q)
{
	if ((q.level == 0) && (fmod(q.index, 1) == 0)) // Exponent is an integer
	{
		SLI Result = exp(log(abs(p)) * q);
		if (p.signi || (p.index < 0))
			return bool(fmod(q.Double(), 2)) ? -Result : Result;
		else
			return Result;
	}

	if ((p.signi==1) || (p.index<0))//negative base
		return SLI(7, 1, 0, 1); // -Inf
	
	// exp(log(*this) * s)
	if ((p.level == 0) && (q.level==0)) // pure FLP
	{
		double Result = pow(p.index, q.index);
		double absResult = ABS(Result);
		if ((absResult > REALMIN) && (absResult < REALMAX)) // not overflow or underflow
			return SLI(Result);
	}
	return exp(log(p) * q);
}

SLI sqrt(const SLI& s)
{
	if (s.level == 0)
		return SLI(sqrt(s.index));
	
	if (s.signi==1)//negative base
		return SLI(7, 1, 0, 1); // -Inf
	
	return exp(log(s) * 0.5);
}

//////////////////////////////////////////////////////////////////////
// Convert a SLI object to a FLP number
//////////////////////////////////////////////////////////////////////
double SLI2FLP(SLI s)
{
	return s.Double();
}

//////////////////////////////////////////////////////////////////////
// Trig functions
//////////////////////////////////////////////////////////////////////
SLI sin(const SLI& s)
{
	double x = s.Double();
	if (ABS(x) < REALMAX)
		return sin(x);
	else if (ABS(x) > 1)
		return 0;
	else
		return s;
}

double cos(const SLI& s)
{
	double x = s.Double();
	if (ABS(x) < REALMAX)
		return cos(x);
	else
		return 1;
}

SLI tan(const SLI& s)
{
	double x = s.Double();
	if (ABS(x) < REALMAX)
		return tan(x);
	else if (ABS(x) > 1)
		return 0;
	else
		return s;
}

//////////////////////////////////////////////////////////////////////
// Inverse trig functions
//////////////////////////////////////////////////////////////////////
double asin(const SLI& s)
{
	return asin(s.Double());
}

double acos(const SLI& s)
{
	return acos(s.Double());
}

double atan(const SLI& s)
{
	return atan(s.Double());
}

//////////////////////////////////////////////////////////////////////
// Hyperbolic trig functions
//////////////////////////////////////////////////////////////////////
SLI sinh(const SLI& s)
{
	double x = s.Double();
	if (ABS(x) <= HYPERMAX)
		return sinh(x);
	else
		return (exp(s) - exp(-s)) * 0.5;
}

SLI cosh(const SLI& s)
{
	double x = s.Double();
	if (ABS(x) <= HYPERMAX)
		return cosh(x);
	else
		return (exp(s) + exp(-s)) * 0.5;
}

double tanh(const SLI& s)
{
	double x = s.Double();
	if (ABS(x) <= HYPERMAX)
		return tanh(x);
	else
		return (exp(s) - exp(-s)) / (exp(s) + exp(-s));
}

//////////////////////////////////////////////////////////////////////
// Quick sort
//////////////////////////////////////////////////////////////////////
void SLIsort(SLI** first, SLI** last)
{
	SLI** left = first;
	SLI** right = last;
	SLI** pivot = left++;
	SLI* temp;
	SLI p, q;

	p = **pivot;
	double pivotAbs = p.index;
	if (pivotAbs < 0)
		pivotAbs = - pivotAbs;

	while (left != right)
	{
		q = **left;
		if ((q.signl > p.signl) ||
		  ((q.signl == p.signl) && (((!p.signl) && ((q.level < p.level) || ((q.level == p.level) && (ABS(q.index) < pivotAbs))) ||
			                        ((p.signl) && ((q.level > p.level) || ((q.level == p.level) && (q.index > p.index))))))))
			left++;
		else
		{
			q = **right;
			while ((left != right) && ((p.signl < q.signl) ||
									 ((p.signl == q.signl) && (((!p.signl) && ((p.level < q.level) || ((p.level == q.level) && (pivotAbs <= ABS(q.index)))) ||
															   ((p.signl) && ((p.level > q.level) || ((p.level == q.level) && (p.index >= q.index)))))))))
				q = **(--right);
			temp = *left;
			*left = *right;
			*right = temp;
		}
	}

	q = **left;
	if ((p.signl < q.signl) ||
	  ((p.signl == q.signl) && (((!p.signl) && ((p.level < q.level) || ((p.level == q.level) && (pivotAbs < ABS(q.index)))) ||
							    ((p.signl) && ((p.level > q.level) || ((p.level == q.level) && (p.index > q.index))))))))
		left--;
	temp = *left;
	*left = *first;
	*first = temp;

	if (first != left)
		SLIsort(first, left);
	if (right != last)
		SLIsort(right, last);
}

//////////////////////////////////////////////////////////////////////
// Summation of N arguments
// Internal operation for a sorted list of SLI numbers
//////////////////////////////////////////////////////////////////////
SLI sum(unsigned N, SLI **s)
{
	register int i, j;
	register int M = N - 1;
	SLI x = *s[M];
	SLI y;
	double a[4], b, c, d, h;
	
	// check if there are more than one number have the largest magnitude
	register unsigned Count = 1;
	i = M - 1;
	y = *s[i];
	while ((i >= 0) && (x.index == y.index) && (x.level == y.level) && (x.signl == y.signl))
	{
		if (x.signi == y.signi)
			Count++;
		else
			Count--;
		i--;
		if (i >= 0)
			y = *s[i];
	}

	i++; // now *s[i] is the first such number
	if (i < M) // there are!
		if (Count == 0) // they cancel each other
			return sum(i, s);
		else // they don't cancel
		{
			x *= Count;
			M = i;
		}

	if (x.level == 0)  // the largest is a FLP number
	{
		double Result = 0;
		for (i = 0; i <= M; i++)
		{
			y = *s[i];
			if (y.level == 0)
				Result += y.index;
			else
				Result += y.Double();
		}
		return SLI(Result);
	}

	if ((x.level > 5) || ((x.level == 5) && (x.index > 0.23125)))  
		// Additions are all trivial after 5.215743659411, but here we allow a little more space
		// so that even 1,000,000 smaller numbers cannot change the largest addend which is bigger than 5.23125
		return x;

	// Generate sequence a
	a[x.level - 2] = exp(-x.index);
	for (j = x.level - 2; j > 0; j--)
		a[j - 1] = exp(-1 / a[j]);

	// Compute c0
	c = 0;
	// compute a0 * b0 ^ r (with sign) for each addend, starting from the largest
	for (i = M - 1; i >= 0; i--)
	{
		y = *s[i];
		//Calculate d
		if (y.level == 0)
		{
			d = ABS(y.index) * exp(-1/a[0]);
			d = y.index < 0 ? -d : d;
		}
		else
		{
			b = exp(-y.index);
			for (j=y.level-2; j>0; j--)
				b = exp(-1/b);
			if (!x.signl)
				if (!y.signl) //large
					d = exp((a[0]-b)/(a[0]*b));
				else //mixed
					d = exp((-a[0]-b)/(a[0]*b));
			else //small
				d = exp((b-a[0])/(a[0]*b));
			d = y.signi ? -d : d;
		}
		c += d;
	}
	c += 1;

	// Generate sequence c
	c = x.signl ? 1 - a[0] * log(c) : 1 + a[0] * log(c);
	j = 0;
	while ((c >= a[j]) && (j < x.level-2))
	       c = 1 + a[++j] * log(c);
	if (c < a[j])
		return SLI(j+1, c/a[j], x.signl, x.signi);
	j++;
	h = x.index + log(c);
	HSEQUENCE
	return SLI(++j, h, x.signl, x.signi);
}

//////////////////////////////////////////////////////////////////////
// Summation of N arguments
// Input number N and a array of SLI numbers
//////////////////////////////////////////////////////////////////////
SLI precisesum(unsigned N, SLI *s)
{
	if (N == 1)
		return *s;

	register unsigned i;
	SLI **ss = new SLI*[N];
	SLI *temp = s;
	for (i = 0; i < N; i++)
		ss[i] = temp++;

	SLIsort(ss, &ss[N-1]);

	SLI Result = sum(N, ss);
	delete[] ss;
	return Result;
}

SLI sum(unsigned N, SLI *s)
{
	if (N == 1)
		return *s;

	register unsigned i;
	SLI **ss = new SLI*[N];
	SLI *temp = s;
	for (i = 0; i < N; i++)
		ss[i] = temp++;

	// find the largest magnitude addend
	SLI MaxAdd = s[0];
	double MaxAbs = ABS(MaxAdd.index);
	unsigned MaxInd = 0;
	unsigned Count = 0;
	unsigned j = N - 1;
	
	while (Count == 0) // do this till find a largest manitude number which will not be canceled by its negative
	{
		N = j + 1;
		for (i = 1; i <= j; i++)
			if ((MaxAdd.signl < s[i].signl) ||
			((MaxAdd.signl == s[i].signl) && (((!MaxAdd.signl) && ((MaxAdd.level < s[i].level) || ((MaxAdd.level == s[i].level) && (MaxAbs < ABS(s[i].index)))) ||
												((MaxAdd.signl) && ((MaxAdd.level > s[i].level) || ((MaxAdd.level == s[i].level) && (MaxAdd.index > s[i].index))))))))
			{
				MaxAdd = s[i];
				MaxInd = i;
				MaxAbs = ABS(MaxAdd.index);
			}

		for (i = MaxInd; i <= j; i++)
			if ((MaxAdd.index == s[i].index) && (MaxAdd.level == s[i].level) && (MaxAdd.signl == s[i].signl))
			{
				if (MaxAdd.signi == s[i].signi)
					Count++;
				else
					Count--;
				ss[j] = &s[i];
				ss[i] = &s[j];
				j--;
			}
	}

	SLI Result = sum(N, ss);
	delete[] ss;
	return Result;
}

//////////////////////////////////////////////////////////////////////
// Summation of N arguments
// Input number N and N arguments separately
//////////////////////////////////////////////////////////////////////
SLI precisesum(unsigned N, SLI start, ...)
{
	if (N == 1)
		return start;

	register unsigned i;
	SLI *pArg = &start;
	SLI **ss = new SLI*[N];
	for (i = 0; i < N; i++)
		ss[i] = pArg++;

	SLIsort(ss, &ss[N-1]);

	SLI Result = sum(N, ss);
	delete[] ss;
	return Result;
}

SLI sum(unsigned N, SLI start, ...)
{
	if (N == 1)
		return start;

	register unsigned i;
	SLI *pArg = &start;
	SLI *s = new SLI[N];
	for (i = 0; i < N; i++)
		s[i] = *pArg++;

	SLI Result = sum(N, s);
	delete[] s;
	return Result;
}

SLI norm(unsigned N, SLI *s)
{
	register unsigned i;
	SLI *x = new SLI[N];

	for (i = 0; i < N; i++)
		x[i] = abs(s[i] * s[i]);
	SLI Result = sqrt(sum(N, x));
	delete[] x;
	return Result;
}

SLI norm(unsigned N, SLI *s, SLI p)
{
	register unsigned i;
	SLI *x = new SLI[N];

	for (i = 0; i < N; i++)
		x[i] = pow(abs(s[i]), p);
	SLI Result = pow(sum(N, x), p.Rec());
	delete[] x;
	return Result;
}

SLI dotproduct(unsigned N, SLI *a, SLI *b)
{
	register unsigned i;
	SLI *x = new SLI[N];

	for (i = 0; i < N; i++)
		x[i] = a[i] * b[i];
	SLI Result = sum(N, x);
	delete[] x;
	return Result;
}

//the output format flag used in ostream
//1 for (level, index); 0 for extended decimal FLP
static bool SLIoutput = 1;

ostream& sliout(ostream& os)
{
	SLIoutput = 1;
	return os;
}

ostream& flpout(ostream& os)
{
	SLIoutput = 0;
	return os;
}

ostream& operator <<(ostream& os, const SLI& s)
{
	if (SLIoutput == 1)
		display(s);
	else
		dispfloat(s);
	return os;
}

std::ostream& operator <<(std::ostream& os, const std::complex<SLI>& c)
{
	return os << "(" << c.real() << ", " << c.imag() << ")";
}
